#include <math.h>
#include <iostream>
#include <sstream>
#include <fstream>

#define NX 32
#define NT 16
#define NAME "CON"
#define PATHIN  "../data/"
#define PATHOUT "../dataOut/"
#define PATHIM  "../img/"
#define NCOOLINGSTEPS 30
#define NFILES 99

using namespace std;

// This file contains four classes: coordinate, matrix, MyLattice and filename.

// coordinate class implements operations with points on a 4-dimentional lattice
// 1 corresponds to x, 2 to y, 3 to z, 4 to t.

// UPDATE: according to Kotow, t becomes 1, and so on.

class filename
{
      public: 
      int number; // just a number from a filename
      string name;
      
      filename(int n);
      string confname();
      string dataname();
      string imagename();
            
      private:
      // current code creates just a number 1, 18, 234.
      string mytostring(int i);
};

class coordinate
{
     public:
     int x, y, z, t;

     coordinate()
     {
     };

     coordinate(int X, int Y, int Z, int T)
     {
          x = X; y = Y; z = Z; t = T;
          CycleFix();
     };
     
     coordinate forw(int dir)
     {
          coordinate cd(x,y,z,t);
          switch(dir)
          {
               case 2: cd.x++; CycleFix(cd.x, NX); break;
               case 3: cd.y++; CycleFix(cd.y, NX); break;
               case 4: cd.z++; CycleFix(cd.z, NX); break;
               case 1: cd.t++; CycleFix(cd.t, NT); break;    
          }
          return cd;
     }
     
     coordinate back(int dir)
     {
          coordinate cd(x,y,z,t);
          switch(dir)
          {
               case 2: cd.x--; CycleFix(cd.x, NX); break;
               case 3: cd.y--; CycleFix(cd.y, NX); break;
               case 4: cd.z--; CycleFix(cd.z, NX); break;
               case 1: cd.t--; CycleFix(cd.t, NT); break;    
          }
          return cd;
     }
     
     coordinate move(int dir)
     {
          coordinate cd(x,y,z,t);
          switch(dir)
          {
               case 2:  cd.x++; CycleFix(cd.x, NX); break;
               case 3:  cd.y++; CycleFix(cd.y, NX); break;
               case 4:  cd.z++; CycleFix(cd.z, NX); break;
               case 1:  cd.t++; CycleFix(cd.t, NT); break;    
               case -2: cd.x--; CycleFix(cd.x, NX); break;
               case -3: cd.y--; CycleFix(cd.y, NX); break;
               case -4: cd.z--; CycleFix(cd.z, NX); break;
               case -1: cd.t--; CycleFix(cd.t, NT); break;    

          }
          return cd;
     }
     
     
     void cprint()
     {
          cout << x << "\t" <<  y << "\t" << z << "\t" << t << "\n";
     }
          
     private:  // Two procedures below are responsible for periodic boundary conditions
     void CycleFix()
     {
          while (x < 0) { x = x + NX;};
          while (y < 0) { y = y + NX;};
          while (z < 0) { z = z + NX;};
          while (t < 0) { t = t + NT;};
          
          x = x % NX;
          y = y % NX;
          z = z % NX;
          t = t % NT;
          
          if(x==0) x = NX;
          if(y==0) y = NX;
          if(z==0) z = NX;
          if(t==0) t = NT;
     };
             
             
     void CycleFix(int &a, int L)
     {
          if(a == L+1) a = 1;
          else if(a == 0) a = L;
     }     
};

class matrix   // The class implements operations with SU(2) matrices in the Pauli matrices representation
{
   public:
     float a0,a1,a2,a3;
     
     matrix()
     {
     };
     
     matrix(float b0, float b1, float b2, float b3)
     {
          a0=b0; a1=b1; a2=b2; a3=b3;
     }
     matrix(const matrix& M)
     {
          a0=M.a0; a1=M.a1; a2=M.a2; a3=M.a3;
     }

     float tr() // Trace
     {
          return 2*a0;
     };
     
     float abssq()  // Determinant
     {
          return a0*a0 + a1*a1 + a2*a2 + a3*a3;
     }

     matrix inv() //Inverse matrix
     { 
          float b0, b1, b2, b3;
          float rho=a0*a0+a1*a1+a2*a2+a3*a3; 
          b0=a0/rho; b1=-a1/rho; b2=-a2/rho; b3=-a3/rho;
          matrix Mat(b0, b1, b2, b3);
          return Mat;
     };  
     
     matrix im()
     {
          float b0, b1, b2, b3;
          b0=0; b1=a1; b2=a2; b3=a3;
          matrix Mat(b0, b1, b2, b3);
          return Mat;            
     };
                  
     void mprint()
     {
          printf("%f %f %f %f\n", a0, a1, a2, a3);
     };

};

     matrix operator+ (matrix A, matrix B) {
         return matrix(A.a0+B.a0,A.a1+B.a1,A.a2+B.a2,A.a3+B.a3); };
     matrix operator- (matrix A, matrix B) {
         return matrix(A.a0-B.a0,A.a1-B.a1,A.a2-B.a2,A.a3-B.a3); };
     matrix operator* (float x, matrix A) {
         return matrix(A.a0*x,A.a1*x,A.a2*x,A.a3*x); };
     matrix operator* (matrix A, float x) {
         return matrix(A.a0*x,A.a1*x,A.a2*x,A.a3*x); };
     matrix operator* (matrix A, matrix B) {
         return matrix(A.a0*B.a0-A.a1*B.a1-A.a2*B.a2-A.a3*B.a3,
               A.a0*B.a1+B.a0*A.a1+A.a3*B.a2-A.a2*B.a3,
               A.a0*B.a2+B.a0*A.a2+A.a1*B.a3-A.a3*B.a1,
               A.a0*B.a3+B.a0*A.a3+A.a2*B.a1-A.a1*B.a2); };
//==============================================================================



class MyLattice
{
      
public:
      float *Data;                               // The 4-dimensional lattice is stored in 1-dimensional
                                                 // array. For the correspondence between lattice coordinates
                                                 // and array number, see function getN(coordinate)
      const static long int N=4*4*NX*NX*NX*NT;
      const static long int NPlaq = 6*NX*NX*NX*NT;
            
      MyLattice()
      {
           Data = new float[N];
           PlaquetteStore = new matrix[NPlaq];
           FmunuStore = new matrix[NPlaq];
      };
      
      ~MyLattice()
      {
           printf("deleting MyLattice");
           delete Data;
           delete PlaquetteStore;
      };
      
	int loadData(FILE* stream)
     {
          int i;
          for(i = 0; i < N; i++)
			if(fscanf(stream, "%f", Data + i)!=1)
			{
				return 1;
			}
		return 0;
     };
      
	void LoadData(filename fname)
     {
		FILE *stream;
		if ( (stream = fopen(fname.confname().c_str(), "r+b")) == NULL)
		{
			cout << "cannot open the file " << fname.confname() << "\n";
			exit(2);
		}
		if(loadData(stream)==1)
		{	
			cout << "cannot read from " << fname.confname() << "\n";
			exit(1);
		}
		fclose(stream);
	};
      
      matrix GetLink(coordinate crd, int dir)
      {
           if(dir > 0)
           {       
              int n = getN(crd) + (dir - 1)*4;
              matrix Mat(Data[n], Data[n+1], Data[n+2], Data[n+3]);
              return Mat;
           }
           else if(dir < 0)
           {
               return GetLink(crd.back(-dir), -dir).inv();
           }
      };
      
      void WriteLink(matrix Lk, coordinate crd, int dir)        
      {
           int n = getN(crd) + (dir - 1)*4;
           Data[n  ] = Lk.a0;
           Data[n+1] = Lk.a1;
           Data[n+2] = Lk.a2;         
           Data[n+3] = Lk.a3;
      };
      
      matrix Plaquette(coordinate crd, int dir1, int dir2) // Plaquette matrix 
      {
           matrix L1 = GetLink(crd, dir1);
           matrix L2 = GetLink(crd.move(dir1), dir2);
           matrix L3 = GetLink(crd.move(dir2), dir1);
           matrix L4 = GetLink(crd, dir2);
           
           return L1*L2*(L3.inv())*(L4.inv());
      };
      
      matrix PlaquetteStored(coordinate crd, int dir1, int dir2)
      {
           if(dir1 == dir2)
           {
                matrix zero(0,0,0,0);
                return zero;  
           }
           else if(dir1 < dir2)
                return PlaquetteStore[getPlaqN(crd, dir1, dir2)];
           else if(dir1 > dir2)
                return PlaquetteStore[getPlaqN(crd, dir2, dir1)].inv();
      }
      
      matrix Fmunu(coordinate crd, int dir1, int dir2)
      {
           if(dir1 == dir2)
           {
               matrix L(0,0,0,0);
               return L;
           }
           
           matrix L1 = Plaquette(crd, dir1, dir2);
           matrix L2 = Plaquette(crd, dir2, -dir1);
           matrix L3 = Plaquette(crd, -dir1, -dir2);
           matrix L4 = Plaquette(crd, -dir2, dir1);
           
           
           return (L1 + L2 + L3 + L4).im()*0.25;
      };
      
      matrix FmunuStored(coordinate crd, int dir1, int dir2)
      {
           if(dir1 == dir2)
           {
                matrix zero(0,0,0,0);
                return zero;  
           }
           else if(dir1 < dir2)
                return FmunuStore[getPlaqN(crd, dir1, dir2)];
           else if(dir1 > dir2)
                return -1*FmunuStore[getPlaqN(crd, dir2, dir1)];
      }
      
      float Tmunu(coordinate crd, int dir1, int dir2)
      {
           float T = 0;
           if(dir1==dir2)
           {
                //It's time to write this part of code 
						
           }
           
           int i;
           for(i = 1; i <=4; i++)
           {
                 if(i != dir1 && i!= dir2)
                 {
                      T+= -(FmunuStored(crd, dir1, i)*FmunuStored(crd, dir2, i)).tr();
                 }
           }
           return T;
      }
      
      float action(coordinate crd, int dir1, int dir2)  // action corresponding to a single plaquette
      {
           return 1 - 0.5*PlaquetteStored(crd, dir1, dir2).tr();
      }
      
      float action(coordinate crd, int dir)             // action corresponding to a single link variable
      {
           float act = 0;
           int dr;
           for(dr = 1; dr <= 4; dr++)
           {
                if(dr == dir) continue;
                act += action(crd, dir, dr);
                act += action(crd.back(dr), dir, dr);
           }
           return act;
      }
      
      
      float Action()                                    // The whole action of the lattice
      {
           float Act = 0;
           int x=0,y=0,z=0,t=0, dir1, dir2;
           coordinate *crd;

                    
           for(x = 0; x < NX; x++)
           for(y = 0; y < NX; y++)
           for(z = 0; z < NX; z++)
           for(t = 0; t < NT; t++)
           {
               crd = new coordinate(x,y,z,t);

               for(dir2 = 1; dir2 <=    4; dir2++)
               for(dir1 = 1; dir1 < dir2; dir1++)
               {
                    Act += action(*crd, dir1, dir2);
               }
               
               delete crd;
           }
                      
           return Act;
      } 
      
// While conducting the cooling procedure, each link is replaced by the result value of the above function

      matrix CoolingStepLink(coordinate crd, int mu)         
      {
           int nu;
           matrix res(0,0,0,0);
           for(nu = 1; nu <= 4; nu++)   
           {
                 if(nu != mu)
                 res = res + 
                       GetLink(crd,          nu)        *
                       GetLink(crd.forw(nu), mu)        *
                       (GetLink(crd.forw(mu), nu).inv())+
                       (GetLink(crd.back(nu), nu).inv())*
                       GetLink(crd.back(nu), mu)*
                       GetLink((crd.back(nu)).forw(mu), nu);
           // The expression above is exactly formula (2) from the arXiv:1401.2441v2 article
           }
                      
           return 1/sqrt(res.abssq())*res;
           
          // cout << 6 - 0.5*res.abssq()*( GetLink(crd, mu) * res.inv() ).tr() << "\n"; 
          // Test for correctness: the value above should be equal to action(crd, mu)  
          // cout << 6 - sqrt(res.abssq()) << "\n";
          // Test for correctness: the value above should be equal to action(crd, mu) after making the cooling elementary step
      };
      
      // Cooling procedure single step as described in arXiv:1401.2441v2
      void CoolingStep()
      {
           int x,y,z,t, dir;
           coordinate *crd;           
           
           
           for(t = 0; t < NT; t++)
           for(x = 0; x < NX; x++)                      
           for(y = 0; y < NX; y++)
           for(z = 0; z < NX; z++)
           {
               crd = new coordinate(x,y,z,t);                  

               for(dir = 1; dir <= 4; dir++)
               {
                    WriteLink(CoolingStepLink(*crd, dir), *crd, dir);
               }                              
               delete crd;
           }
           
//           StorePlaquettes();
      }

      void CoolingStepInv()
      {
           int x,y,z,t, dir;
           coordinate *crd;           
           
           
           for(t = 0; t < NT; t++)
           for(x = NX-1; x >= 0; x--)                      
           for(y = 0; y < NX; y++)
           for(z = 0; z < NX; z++)
           {
               crd = new coordinate(x,y,z,t);                  

               for(dir = 1; dir <= 4; dir++)
               {
                    WriteLink(CoolingStepLink(*crd, dir), *crd, dir);
               }                              
               delete crd;
           }
           
//           StorePlaquettes();
      }
      
      void CoolingStepInv2()
      {
           int x,y,z,t, dir;
           coordinate *crd;           
           
           
           for(t = 0; t < NT; t++)
           for(x = NX-1; x >= 0; x--)                      
           for(y = 0; y < NX; y++)
           for(z = 0; z < NX; z++)
           {
               crd = new coordinate(x,y,z,t);                  

               for(dir = 3; dir >= 1; dir--)
               {
                    WriteLink(CoolingStepLink(*crd, dir), *crd, dir);
               }                              
               WriteLink(CoolingStepLink(*crd, 4), *crd, 4);               
               delete crd;
           }
           
//           StorePlaquettes();
      }

      void CoolingStepInv3()
      {
           int x,y,z,t, dir;
           coordinate *crd;           
           
           
           for(t = 0; t < NT; t++)
           for(x = NX-1; x >= 0; x--)
           for(y = 0; y < NX; y++)
           for(z = 0; z < NX; z++)
           {
               crd = new coordinate(x,y,z,t);                  

               for(dir = 4; dir >= 1; dir--)
               {
                    WriteLink(CoolingStepLink(*crd, dir), *crd, dir);
               }                              
               delete crd;
           }
           
//           StorePlaquettes();
      }

      
      // local      
      float topologicalCharge(coordinate crd)
      {
            return  (PlaquetteStored(crd, 1,2)*PlaquetteStored(crd, 3,4)).tr()
                  - (PlaquetteStored(crd, 1,2)*PlaquetteStored(crd, 4,3)).tr()
                  - (PlaquetteStored(crd, 1,3)*PlaquetteStored(crd, 2,4)).tr()
                  + (PlaquetteStored(crd, 1,3)*PlaquetteStored(crd, 4,2)).tr()
                  + (PlaquetteStored(crd, 1,4)*PlaquetteStored(crd, 2,3)).tr()
                  - (PlaquetteStored(crd, 1,4)*PlaquetteStored(crd, 3,2)).tr();
      };
      // global
      float TopologicalCharge()
      {
           float Charge = 0;
           int x=0,y=0,z=0,t=0;
           coordinate *crd;

                    
           for(x = 0; x < NX; x++)
           for(y = 0; y < NX; y++)
           for(z = 0; z < NX; z++)
           for(t = 0; t < NT; t++)
           {
               crd = new coordinate(x,y,z,t);
               Charge += topologicalCharge(*crd);
               delete crd;
           }
           
           return Charge;
      };


      void StorePlaquettes()
      {
           int x,y,z,t, dir1, dir2;
           coordinate *crd;
           
           for(x = 0; x < NX; x++)
           for(y = 0; y < NX; y++)
           for(z = 0; z < NX; z++)
           for(t = 0; t < NT; t++)
           {
               crd = new coordinate(x,y,z,t);

               for(dir2 = 1; dir2 <= 4; dir2++)
               for(dir1 = 1; dir1 < dir2; dir1++)
               {
                    PlaquetteStore[getPlaqN(*crd, dir1, dir2)] = Plaquette(*crd, dir1, dir2);
               }
               
               delete crd;
           }                      
      }

      void StoreFmunu()
      {
           int x,y,z,t, dir1, dir2;
           coordinate *crd;
           
           for(x = 0; x < NX; x++)
           for(y = 0; y < NX; y++)
           for(z = 0; z < NX; z++)
           for(t = 0; t < NT; t++)
           {
               crd = new coordinate(x,y,z,t);

               for(dir2 = 1; dir2 <= 4; dir2++)
               for(dir1 = 1; dir1 < dir2; dir1++)
               {
                    FmunuStore[getPlaqN(*crd, dir1, dir2)] = Fmunu(*crd, dir1, dir2);
               }
               
               delete crd;
           }
           cout << "Fmunu stored\n";    
      }


      
private:
      matrix *PlaquetteStore;
      matrix *FmunuStore;
      
      long int getN(coordinate crd)
      {
           return 16*( (crd.t-1) + (crd.x-1)*NT + (crd.y-1)*NT*NX + (crd.z-1)*NT*NX*NX);
      };

      long int getPlaqN(coordinate crd, int dir1, int dir2)
      {
           int s;
           if(dir1 == 1 && dir2 == 2)       s = 0;
           else if(dir1 == 1 && dir2 == 3)  s = 1;
           else if(dir1 == 1 && dir2 == 4)  s = 2;
           else if(dir1 == 2 && dir2 == 3)  s = 3;
           else if(dir1 == 2 && dir2 == 4)  s = 4;
           else if(dir1 == 3 && dir2 == 4)  s = 5;
           else exit(1);
           return 6*(NX*NX*NX*(crd.t-1) + NX*NX*(crd.z-1) + NX*(crd.y-1) + (crd.x-1)) + s;
      };
      
};


filename::filename(int n)
{
     number = n;
     name = NAME + mytostring(n);
}

string filename::confname()
{
     return (PATHIN + name + ".LAT").c_str();
}

string filename::dataname()
{
     return (PATHOUT + name).c_str();
}
      
string filename::imagename()
{
     return (PATHIM + name + ".jpg").c_str();
}
            
string filename::mytostring(int i)
{
     stringstream ss;
     ss << i;
     string num = ss.str();
             
     return num;
}

